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Continuum hydrodynamic models of active liquid crystals have been used to describe dynamic self- 
organising systems such as bacterial swarms and cytoskeletal gels. A key prediction of such models 
is the existence of self-stabilising kink states that spontaneously generate fluid flow in quasi-one 
dimensional channels [l|]. Using simple stability arguments and numerical calculations we extend 
previous studies to give a complete characterisation of the phase space for both contractile and 
extensile particles (ie pullers and pushers) moving in a narrow channel as a function of their flow 
alignment properties and initial orientation. This gives a framework for unifying many of the results 
in the literature. We describe the response of the kink states to an imposed shear, and investigate 
how allowing the system to be polar modifies its dynamical behaviour. 

PACS numbers: 87.10.-|-e, 83.80.Xz, 47.60.-i 



Suspensions of flagellate bacteria swim in quasi- 
turbulent swirls that suggest long-range collective order- 
ing f5]. Cytoskeletal filaments and motor-proteins spon- 
taneously form stars and spirals without any external 
guide Q. These are just two examples of a class of self- 
organizing system that poses a novel challenge to physi- 
cists, as the ordering is intrinsically non-equilibrium in 
nature. The pattern formation is a dynamic phenomenon 
that relies on a continuous expenditure of energy by the 
individual particles as they actively generate forces on 
each other and/or the surrounding medium 

Convergent strands of research focused on modeling 
bacterial swarms 0, and on cytoskeletal dynam- 
ics [13, 11 , Hi have suggested that the collective be- 
haviour seen in both types of system can be described in 
the continuum limit by the same phenomenological, hy- 
drodynamic model. This builds on the equations of liquid 
crystal hydrodynamics, which capture the inherent di- 
rectionality of the particles (e.g. bacteria, microtubules), 
and adds extra non-equilibrium terms to account for the 
activity. One of the most striking predictions of the ac- 
tive liquid crystal model has been that, under certain 
conditions, the uniform aligned state is unstable |5i]; fur- 
thermore, when confined to a quasi-one dimensional slab 
geometry by flat interfaces (fig.H]) this instability can lead 
to the onset of a spontaneous steady flow |ll,[l3| • Previous 
studies of the spontaneous flow transition in the active 
liquid crystal model using both analytic f]\ and numerical 
[i4.,i5, 16] approaches have focused on illuminating only 
certain neighbourhoods of the complete phase space. To 
help map the phenomenological coefficients in the con- 
tinuum theory onto microscopic parameters, and to aid 
in identifying their values for a given physical system, 
it is helpful to have a more complete picture of the ef- 
fect of varying key coefficients. Therefore our aim in the 
Letter is to extend and unify previous studies by carry- 
ing out a complete numerical exploration of the phase 
space formed by three of the most important physical 
characteristics of the system: the type and strength of 
the flow field generated by the individual particles, their 
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FIG. 1: Quasi-one dimensional slab geometry, with solid, no- 
slip walls lying parallel to the i-axis ensuring that only the 
x-component of the fluid velocity is non-zero. The orientation 
angle of the director field is measured such that S = is 
parallel to the walls, and 6 = n /2 is perpendicular. 



fiow alignment behaviour, and the initial orientation of 
the director field. We distinguish six spontaneous flow 
states, describe their response to an applied shear, and 
explain how they are modified in a system with polar 
symmetry. 

We assume that, apart from their activity, the assem- 
bly of particles can be well approximated as nematic liq- 
uid crystal — that is, there is no difference between the 
two ends of an individual particle. The effect of relaxing 
this assumption by allowing the particles to be polar is 
discussed at the end of this letter. The hydrodynamics of 
nematic liquid crystals is well developed [13] ■ We adopt 
the Erickson-Leslie-Parodi approach [1], in which the ne- 
matic order parameter is a fixed-magnitude unit vector 
field fi which evolves according to 
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where u is the velocity field of the solvent; A is the flow 
alignment parameter, more on which below; T is a rota- 
tional viscosity; h is the molecular field given by 
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FIG. 2: The dipolar flow flelds (curved arrows) generated by 
extensile (left) and contractile (right) particles. The straight 
vertical arrows represent the director fleld, which is along the 
rod axis for rod-like (A > 0) particles and perpendicular to 
the disk plane for discoidal (A < 0) particles. 

(assuming only a single elastic constant K); and 

= ^{diUj + djUi) ; flij = ^{diUj - djUi). (3) 

The first two terms on the right-hand side of eq. ([T]) de- 
scribe alignment (or tumbling) of the director field by 
local shear flow. The third term accounts for the ten- 
dency of the ordered nematic to resist distortions, and 
arises ultimately from excluded volume interactions be- 
tween individual particles. 

The flow field u obeys the Navier-Stokes equation 

p{dt + u ■ V)ui = dj{aij) + rjdj{diUj + d-jUi) (4) 

where p is the fluid density and t] is the viscosity. The 
passive part of the stress tensor o-y is 

crfj = -p6rj - ^[n^hJ + Ujh,] + ^[riihj - njhi] (5) 

where p is the bulk pressure. 

Eqs. (1) — (5) are generalised to describe active sys- 
tems by introducing one or more terms cr° to the stress 
tensor, so that a = a^ + a"'. These active stresses account 
for the forces generated by the individual particles, as a 
function of the local director field, and cannot be derived 
from a free energy; however, their form can be arrived at 
from considering the symmetry of the flow field generated 
by the particles. It is commonly assumed that this flow 
field can be described to leading order as dipolar, and 
that higher order contributions can be neglected, when 
considering far-fleld interactions and long time scales Q . 
The appropriate active stress tensor for a suspension of 
force dipoles is 

where C is the activity coefficient, the meaning of which 
is discussed below. One can also add active terms to the 
director equation of motion, eq. ([T]), to account for self- 
alignment effects fl[, but as such terms do not play an 



important role in the spontaneous flow transition they are 
neglected here. For the same reason, we neglect active 
contributions to the isotropic pressure p. 

The two key parameters that characterise the individ- 
ual particles are the activity parameter ^ and the flow 
alignment parameter A. Both the magnitude and sign 
of each have an easily understood physical significance. 
The sign of ^ determines whether the dipolar flow field 
generated by the particles is extensile (C > 0) or contrac- 
tile (C < 0), as illustrated in fig. El In the swimmer lit- 
erature, an alternative nomenclature is sometimes used, 
with extensile swimmers described as pushers and con- 
tractile swimmers as pullers. The magnitude of C relates 
to the strength of the forcing, or equivalently the rate at 
which the particles expend energy. 

The flow alignment parameter A determines how the 
director field responds to a shear flow. A shear flow can 
be decomposed into an extensional and a rotational com- 
ponent, and the magnitude of A determines their rela- 
tive influence. For |A| < 1, the rotational part of the 
flow always dominates regardless of the director orienta- 
tion, and thus the director will continuously rotate under 
shear. This is the flow tumbling regime. For jA| > 1, 
however, the director will tend to align at a unique angle 
to the flow direction, 9l = ^ cos~^ at which the ex- 
tensional and rotational parts of the shear flow balance. 
This is the flow aligning regime. In general, molecule- 
sized nematic liquid crystals are flow aligning and larger 
objects are flow tumbling, since for the extensional flow 
to dominate there must be enough Brownian rotational 
diffusion to wash out the rotational part of the flow. The 
sign of A is also important, and relates to the shape of the 
individual particles. Rod-like objects have A > 0; discoid 
objects have A < 0; and A = corresponds to spheres. 

The original analytic prediction of the spontaneous 
flow transition by Voiturez et al. [ij considered the con- 
tractile (C < 0), discoidal, flow aligning (A < —1) portion 
of phase space. Previous numerical studies by Maren- 
duzzo et al. considered both extensile and contractile 
systems, for rod-like particles (A > 0). We numerically 
explore the complete (C, A) phase space, and vary the 
orientation 0o of the initial uniformly aligned state, to 
access a number of novel spontaneous flow states. 

Eqs. (l)-(6) were solved, using a hybrid Lattice Boltz- 
mann approach analogous to that described by Maren- 
duzzo et al. [isj , in the one dimensional geometry of Fig. 
[H We follow Gates et al. in using free boundary con- 
ditions for the director fleld {dyUi = at the walls); the 
effect of fixed boundary conditions will be mentioned be- 
low. The following parameter values were used through- 
out: slab width L = 100, K = 0.04, rj = 1.27, and 
r — 0.2 (all in simulation units). The director field was 
initially uniform and aligned at an angle 9q to the walls. 
To break the symmetry we perturbed the uniform state 
by generating a small kink at the centre of the slab by ro- 
tating the director in the left half of the domain clockwise 
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FIG. 3: Stable states (in activity ( versus flow alignment parameter A) for a one-dimensional slab of active nematic. Initial 
director orientation (a) parallel to the walls 9o — 0. (b) perpendicular to the walls So = f- In the shaded (unshaded) regions 
the initial state is unstable (stable) to small perturbations. The phase boundaries are given by eq. ([7]). The diagrams show the 
configuration of the director field n for the steady spontaneous flow state that arises in each unstable region, (c) Example flow 
fields Ux{y) corresponding to the spontaneous flow states identified in (a) and (b). Parameters used are: for El/Cl, A = +1.5 
/ -1.5, C = +0.001/ -0.001, 6^0 = / f ; for E2/C2, A = +1.5 / -1.5, ( = +0.001/ -0.001, 6»o = f / 0; and for E3/C3, A = / 0, 
= +0.001/ -0.001, So = f- / 0. Other parameters are given in the text. 
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FIG. 4: Response of an El-type kink to an applied shear, and the corresponding induced fiow profile. The shear is generated 
by moving the right wall with a steady upwards velocity Hwaii- From left to right, Uwaii ~ 0.005, 0.01, 0.02 in simulation units. 
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0.1 radians and that in the right half the same amount 
anti-clockwise. The diagrams in fig. [3] show where the 
system was found to be stable (unshaded) and unstable 
(shaded) against the perturbation, for (a) 6*0 = 0, and 
(b) 6*0 — tt/2. The critical value of the activity Cc cor- 
responding to the phase boundary is determined analyt- 
ically by linearizing the equations of motion around the 
initial state, and checking when a non-uniform solution 
exists, following Voituriez et al. [H 
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The director configurations in fig. [3] are typical numer- 
ical results obtained after the spontaneous flow state is 
reached in each of the unstable regions. There are six in 
total, three each for extensile (labeled El, E2, E3) and 
contractile (CI, C2, C3) systems. In the previous work 
by Marenduzzo et al. [3, [l3|, regions El and C3 were 
explored. The C2 state was predicted analytically in Voi- 
turiez et al. [ij but has not yet been studied numerically. 

The six kink states are, for the geometry we consider, 
related in pairs by a symmetry transformation. This 
symmetry is apparent from an examination of the flow 
profiles generated by each type of spontaneous flow state, 
shown in fig. ^c). Notice that the states naturally fall 
into three pairs (El/Cl, E2/C2 and E3/C3), each with 
the same type of flow profile, but director fields rotated 
by 7r/2 radians with respect to each other. The states in 
each pair have ^ and A values with the same magnitudes, 
but opposite signs. This highlights the symmetry inher- 
ent in the equations of motion: a change in sign of ( is 
equivalent to a change in sign of A together with a, tt/2 
rotation of the director field; an extensile rod-like particle 
interacts with a flow field the same way as a contractile 
discoidal particle rotated by 7r/2. Note that although 
this holds true for the quasi-one dimensional geometry 
of interest, it is not universally true; in particular the 
symmetry does not hold in three dimensions. 

The instabilities leading to each type of kink can be un- 
derstood from a simple pictorial model, extending ideas 
from T9] (who consider the fiow tumbling regime). Fig. 
[5] illustrates kink perturbations to the initial state of the 
director field, = on the left and Oq — tt/2 on the 
right. The discontinuity in n at the centre of the kink 
generates a force on the fluid due to the active term in 
the stress tensor <tP. The flow direction depends on the 
sign of C; in the figure, all kinks are drawn such that the 
resulting flow is upwards. Because of the no-slip condi- 
tion at the walls, this central forcing necessarily generates 
a region of uniform shear flow on either side of the kink, 
illustrated by the fat vertical arrows. The wide arrows 
overlaid on the director field show the effect of the shear 
fiow on the director in each case: depending on the val- 
ues of A and either the extensional fiow dominates, at- 



tempting to align the director at the Leslie angle (double 
headed arrows) or the rotational flow dominates (curved 
arrows) and causes the director to tumble. By reacting 
thus to the shear fiow, the director field either returns to 
the initial aligned state and the induced shear fiow van- 
ishes, or the perturbation grows until one of the stable 
spontaneous fiow states is reached. The same instability 
mechanism applies in two dimensions, except that with- 
out walls to restrict fiow to one axis there will be no 
steady spontaneous flow states. Instead, the instabilities 
drive quasi-turbulent mixing flow 0, [§] . 

The six stable kink states in flg.[3]are the fundamental 
building blocks for more complicated states observed in 
other numerical studies. For instance, if flxed boundary 
conditions are used for the director field as in Maren- 
duzzo et al. [lit , a half-kink is necessary at each wall. At 
higher activities, striped states can be observed, which 
are simply multiple kinks sitting side by side. We do 
not observe striped states in our study, even at higher 
CI values, because our method for perturbing the ini- 
tial state seeds a single kink at the centre. In contrast, 
Marenduzzo et al. perturb by increasing the director ori- 
entation slightly at a single point, which stimulates the 
formation of at least two kinks, one on either side of the 
perturbation. Thus the exact state chosen by the system 
can depend strongly on how the initial uniform state is 
perturbed. It is possible to have the two different types 
of flow aligning kink (El and E2, or CI and C2) coex- 
isting in the same state, but only if the state is specially 
prepared. Note that we restrict the director to the x-y 
plane, which does not allow for twisted states, as were 



sometimes observed for high |C| in 15|. 

Each of the states reacts differently to externally im- 
posed shear. For the tumbling states, E3 and C3, even 
a small amount of imposed shear causes the director to 
rotate continuously, and there is no steady state. Indeed 
these states are only stable in the first place because they 
spontaneously adopt a configuration which generates no 
fiuid fiow, and thus no shear, outside the kink. In con- 
trast, the flow aligning states El/Cl and E2/C2 can sur- 
vive in applied shear up to a certain magnitude. A single 
kink of either type responds to applied shear by moving 
towards the wall that is moving in the same direction as 
the spontaneous flow. This is illustrated in fig. [4] for an 
El kink. At some critical wall velocity, the kink gets too 
close to the wall and becomes unstable; for higher veloci- 
ties, there is no steady kink state and the velocity profile 
is standard linear shear flow. 

Having catalogued this variety of states, it is impor- 
tant to ask which we would predict to occur in specific 
physical systems. Considering first swimming bacteria: 
they may be either extensile or contractile, depending 
on their mode of swimming, but most are certainly rod- 
like (A > 0) and furthermore should be flow tumbling 
(|A| < 1) given their relatively large size. Therefore, with 
reference to fig. [31 one might expect to see states of type 
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FIG. 5: Mechanism driving the (in)stability of an initiaUy uniformly ahgned state. Thin straight arrows show the director 
field, following a kink perturbation acting on a state with initial orientation So = (left) and do = 7r/2 (right). Arrowheads 
pointing outwards (inwards) denote extensile (contractile). Wide vertical arrows show shear flow generated by the kink. Wide 
arrows superimposed on the director field show the influence of the shear flow upon it: either aligning (straight, double-headed 
arrows) or tumbling (curved arrows). 



E3 for suspensions of "pushers" (e.g. E. coli) or C3 for 
"pullers" (e.g. Chlamydomonas) . This is consistent with 
the experiments and analytical predictions of Berke et al. 
[20| showing that pushers tend to align parallel to walls, 
as seen in state E3, while pullers align perpendicularly as 
in C3. Note however that the existence of spontaneous 
flow in these states relies on a non-zero elastic constant 
K; if K 0, the width of the kink also goes to zero, 
and for the states E3 and C3 the induced flow would 
vanish. Therefore attempts to simulate suspensions of 
microswimmers must incorporate liquid crystalline elas- 
ticity, for example via excluded volume interactions, to 
see steady flow-generating kinks of this type. Previous 
simulations of large numbers of rod-like swimmers in the 
literature 
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|22| | do not incorporate excluded volume 
interactions, and therefore would not be expected to re- 
produce any such steady spontaneous flow states. How- 
ever, our predictions are consistent with the instabilities 
of uniform aligned states observed in these simulations. 
In particular, the enhancement of mixing by pushers (^ex- 
tensile rods) but not by pullers (contractile rods) [1,[2^ 
can be explained with reference to fig. [51 For the value 
A = 1 assumed by Saintillan et al. |22| . an aligned state 
of contractile rods is stable whereas one of extensile rods 
is unstable. 

For solutions of cytoskeletal rods and molecular motors 
there is experimental evidence for the active stress being 
contractile [ij], but it is not clear whether they should 
be treated as flow aligning or flow tumbling. Voituriez 
et al. [l| assumed the former, in which case the particles 
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FIG. 6: The twelve distinct kink types that can occur when 
the particles are polar. Only six of them — those with 
'up/down' alignment — are persistent long-term states, with 
the others moving sideways after formation until they are an- 
nihilated at a wall. 



must be discoidal (A < — 1) for spontaneous flow to occur, 
via states CI or C2. If they are flow tumbling, then they 
could be expected to generate kinks of type C3. In either 
case, our results suggest that cytoskeletal components 
will tend to align roughly perpendicular to the walls, if 
the director is not otherwise pinned at the boundary. 

The story changes significantly if the particles are al- 
lowed to be polar, having a preferred direction along 
which they will move. Polarity is accounted for in the 
model by including more active terms that are forbidden 
by symmetry in the nematic case. At lowest order, these 
comprise an extra active contribution to the stress tensor 
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af° = (3{dinj + djUi), and the replacement of u in the 
director equation of motion (eq. [T|) with u + f3n (for sim- 
phcity we assume that both these terms carry the same 
polar activity coefficient /3, but in general they may be 
different). When /? is allowed to be non-zero, the num- 
ber of distinct kink states is doubled with respect to the 
nematic case (fig. [6]) since each of the original kinks il- 
lustrated in fig. [3] now has two possible alignments. The 
nematic states for which the kink centre is parallel to the 
walls (El, C2 and C3) each have 'up' and 'down' polar 
analogues; the nematic kinks with centres perpendicular 
to the walls (CI, E2 and E3) have 'left' and 'right' ana- 
logues. Numerical studies of the polar kinks show that 
only the 'up/down' kinks remain stationary and persis- 
tent. The 'left/right' kinks, by contrast, move in the 
direction they point until they reach a wall, and then 
vanish, leaving a uniformly aligned state with zero in- 
duced flow. Note that the symmetry between the three 
pairs of states El/Cl, E2/C2 and E3/C3 is broken by the 
introduction of polarity. The remaining stationary kinks 
can be further subdivided into two groups: those that 
point inwards (Elup, G2down and CSdown), for which 
the sharpness of the peak in the induced velocity profile 
increases with /3, and those that point outwards {Eldown, 
C2up and C3up) for which it decreases. If the local con- 
centration of particles is allowed to vary, evolving via 
diffusion and self-advection, the particle concentration is 
enhanced at the centre of the kink and reduced at the 
walls for the inward kinks, but tends to build up at the 
walls for the outward kinks. 

We have unified and extended previous studies 
of spontaneous flow in quasi- ID active nematics by 
numerically exploring the range of possible steady 
states as the activity, flow alignment characteristics and 
orientation are varied. We predict three qualitatively 
distinct types of spontaneous flow state, and examine 
how they respond to externally applied shear. We have 
also studied how the ensemble of stable states is affected 
when the active particles are made polar. Knowledge of 
the range of possible states should help match continuum 
and microscopic models of active materials and guide 
future experiments and simulations attempting to realize 
spontaneous fluid flow in narrow channels. 
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